Triangular lattice quantum dimer model with variable dimer density

Quantum dimer models are known to host topological quantum spin liquid phases, and it has recently become possible to simulate such models with Rydberg atoms trapped in arrays of optical tweezers. Here, we present large-scale quantum Monte Carlo simulation results on an extension of the triangular lattice quantum dimer model with terms in the Hamiltonian annihilating and creating single dimers. We find distinct odd and even \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathbb{Z}}}_{2}$$\end{document}Z2 spin liquids, along with several phases with no topological order: a staggered crystal, a nematic phase, and a trivial symmetric phase with no obvious broken symmetry. We also present dynamic spectra of the phases, and note implications for experiments on Rydberg atoms.

Recent quantum simulation advances have provided remarkable microscopic access to the quantum correlations of a Z 2 quantum spin liquid (QSL) 1,2 . The Z 2 QSL 3,4 is the simplest quantum state in two spatial dimensions with fractionalized excitations and time-reversal symmetry, and has the same anyon content as the toric code 5 . Once we include considerations of lattice and other symmetries, Z 2 QSLs come in different varieties; the distinctions between them are important in understanding the phase diagrams of possible experimental realizations. The coarsest classification subdivides Z 2 QSLs into "odd" and "even" classes, depending upon whether elementary translations anticommute or commute when acting on excitations carrying Z 2 magnetic flux [6][7][8][9] , and results in different translational symmetry fractionalization patterns and spectral signatures in the dynamic response [10][11][12][13][14][15] . More refined classifications have been obtained since [16][17][18][19][20] .
Quantum dimer models (QDMs) 21,22 on nonbipartite lattices have long been known to host Z 2 QSLs. In this work, we investigate an important-but hitherto unexplored-extension of the quantum dimer model on the triangular lattice [23][24][25] . Unlike the more conventionally studied QDMs, here, the density of dimers is allowed to vary by terms in the Hamiltonian which can annihilate and create single dimers on each link of the triangular lattice. Such a dimer-nonconserving term is motivated by connections to models of ultracold atoms trapped in optical tweezers 26,27 , in which each dimer is identified with an atom excited to a Rydberg state by laser pumping [28][29][30] . The observations of ref.
2 are for the case where the atoms are positioned on the links of the kagome lattice; this connects to the quantum dimer model on the kagome lattice 29 . Our study pertains to the triangular-lattice dimer model, which connects to the case where the atoms are placed on the sites of the kagome lattice 24,25,28 ; such a configuration can be readily realized in the experiments, and initial explorations of quantum phases in such a lattice have already been carried out by the team of ref. 2.
With a dimer-nonconserving term present, here we show, the triangular-lattice quantum dimer model displays novel features relevant to the Rydberg-atom experiments. When the nonconserving terms are large, we can obtain a 'trivial' phase with neither topological order nor broken lattice symmetry. More interestingly, the phase diagram of this extended QDM also harbors both odd and even Z 2 liquids. Note that in early discussions of such QSLs in dimer models, the distinction between the liquids was tied to whether the number of dimers on each site was constrained to be odd or even 24,25 . In the present model, the number of dimers on each site fluctuates between odd and even values, namely 1 and 2; nevertheless, the distinction between even and odd QSLs still survives based on the symmetry transformation properties of excitations with magnetic Z 2 flux ("visons"). In the case with a dimer number constraint on each site, there is an anomaly relation requiring that odd (even) dimers produce vison translations that anticommute (commute) 18,19 . However, in the case without a dimer number constraint (or a soft constraint), of interest to us here, microscopic details will determine whether vison translations anticommute or commute, and we will investigate this fate numerically with quantum Monte Carlo simulations.
Finally, our study also obtains several phases which break lattice symmetries, but are topologically trivial. This includes two "staggered" phases 23 , a "columnar" phase 31 , and a "nematic" phase 24,25,32 , and we also discuss their density-wave-ordered counterparts in the context of experiments on Rydberg quantum simulators.

The model
We investigate the following general dimer Hamiltonian, with one or two dimer(s) per site, on the triangular lattice, where the sum on r runs over all plaquettes (rhombi), including the three possible orientations, and l runs over all links. The different terms in this Hamiltonian are as follows. The kinetic term (controlled by t) flips the two dimers on every flippable plaquette, i.e., on each plaquette with two parallel dimers, while the potential term (controlled by the interaction V) describes a repulsion (V > 0) or an attraction (V < 0) between nearest-neighbor dimers. The transversefield term of strength h creates/annihilates a dimer at link l (similar terms also appear in the quantum realization of the classical models of ref. 32), in contrast to the t and V terms, neither of which change the dimer number. Lastly, μ sets the chemical potential for the occupation of a link by a dimer. We further impose a soft constraint requiring that there must be one or two dimer(s) per site. Thus, when μ → ±∞, the model reverts to the conventional hard-constrained quantum dimer model with exactly two or one dimer(s) per site-the phase diagrams of both these QDMs have been extensively studied in the literature 24,25,31,[33][34][35][36][37] . Hereafter, we set t = 1 as the unit of energy for the rest of this paper.
To solve the model in Eq. (1) in an unbiased manner, we employ the recently developed sweeping cluster quantum Monte Carlo algorithm, which can perform efficient sampling in constrained quantum many-body systems [37][38][39][40] . By monitoring the behavior of various physical observables such as dimer correlation functions and structure factors, we map out the detailed phase diagrams, such as, for instance, in Fig. 1. Moreover, in addition to static observables, we also compute the dynamic dimer correlation functions in imaginary time and employ the stochastic analytic continuation method 12,13,37,[41][42][43][44][45][46] to obtain the dynamic dimer spectral functions in real frequencies. Our simulations are performed on the triangular lattice with periodic boundary conditions and system sizes N = 3L 2 for linear dimensions L = 8, 12, 16, 18, 24, while setting the inverse temperature β = L (β = 200) for equal-time (dynamical) simulations.

The phase diagram
Although the phase diagrams in the two limits with exactly 1/3 and 1/6 dimer fillings are well understood, the manner in which they connect to each other in the presence of a nonzero transverse field h and chemical potential μ is an interesting open question. In particular, one may ask what happens between the two kinds of Z 2 QSLs, i.e., whether they are separated by a direct phase transition or an intermediate phase. An important reason this question has remained unaddressed so far is the lack of a suitable algorithm to deal with the soft constraint. As discussed in detail in the section in "Methods", here, we adapt the sweeping cluster Monte Carlo algorithm used for hard-constrained QDMs 38,39 to soft ones and use it to map out the phase diagram of the Hamiltonian in Eq. (1). Figure 1 shows the full phase diagram obtained at h = 0.4, which we focus on in the main text, leaving the discussion of similar phase diagrams with different h to Supplementary Notes 2 and 3 of the Supplementary Information (SI).
The phase diagram exhibits four different symmetry-breaking phases, including the nematic, the columnar, and two staggered phases; the schematic plots of these crystalline phases are shown in the right panels of Fig. 1. Furthermore, we observe two distinct Z 2 QSL phases, which are denoted as "Even QSL" and "Odd QSL" in the figure.
In addition, a trivial disordered-or paramagnetic (PM)-phase exists in the central region in between the two QSLs; note that such a PM phase -columnar) or first-order (such as the QSL--staggered). Right panel: Schematic pictures of the four crystalline phases (nematic, columnar, 1/3 staggered, and 1/ 6 staggered). In the limit of exactly one dimer per site, a ffiffiffiffiffi 12 p × ffiffiffiffiffi 12 p valence bond solid (VBS) phase is known to exist between the odd QSL and the columnar phase. However, it is nearly degenerate with the columnar phase over a large region in our simulations, and we depict this schematically by using a lighter shading for the columnar phase near the odd QSL.
does not arise in the more familiar QDMs where the dimer number per site is exactly constrained. The phase boundaries between these phases are determined by examining various parameter points and paths scanning through the phase diagram, such as the dashed line in Fig. 1.
To characterize this rich variety of phases, we compute the equaltime (τ = 0) dimer-structure factor (see Fig. 2) as where n i is the dimer number operator on bond i and α stands for the three bond orientations, at five representative parameter points corresponding to the five different phases in the phase diagram. Figure 2a-c shows D(k, 0) inside the odd QSL, PM, and even QSL phases, respectively. In the hexagonal Brillouin zone, we observe that there are no peaks associated with long-range order but only broad profiles signifying different short-range dimer correlation patterns in real space. In contrast, Fig. 2d, e presents the dimer-structure factors inside the columnar and nematic phases, respectively. One now clearly sees the Bragg peaks at the M points for the columnar phase (there can be three different orientations of the columnar dimers, corresponding to all the three pairs of M points), and at the Γ point in the nematic phase.

The two Z 2 QSLs
Having established the lack of long-range dimer-dimer correlations in the odd/even Z 2 QSLs and the PM phase, next, we move on to the phase transitions between them. Since all three of these phases are disordered, care needs to be taken in determining their phase boundaries. Our results in this regard are summarized in Fig. 3, which shows the data along a path with a fixed V = 0.9 and varying μ in the phase diagram (dashed line in Fig. 1). First, in Fig. 3a, we illustrate the energy density curves, which appear to be smooth without any obvious turning points along the path as μ is scanned. However, when the transverse field becomes large, we expect that all the links should be polarized along the x axis (if there were no constraints). Since the model in Eq. (1) can be regarded as a spin model with spins on links (occupied/empty links being equivalent to spin up/down), the polarization can be used to describe the level of polarized links (spins), and thus, to probe the PM phase. Indeed, as seen in Fig. 3b, M x helps us to identify a first-order phase transition between the PM phase and the two Z 2 QSLs. In the PM phase, M x becomes large but is still far from the classical saturation value of 1; this is because the soft constraint forbids all links from being fully polarized simultaneously. We can also discover similar first-order phase transitions, at the same parameter points, independently from the dimer filling ρ shown in Fig. 3c. In the even (odd) Z 2 QSL phase, the filling is nearly 1/3 (1/6) while the filling changes continuously in the PM phase.
In addition, a closed string operator 2 , schematically defined as in Fig. 3e, f as 〈string〉 = 〈(−1) #cut dimers 〉 on a rhomboid with odd linear size, can be used to distinguish the two QSLs and the PM phase. As shown in Fig. 3e, f, 〈string〉 should be ±1 in a pure even/odd Z 2 QSL without spinons and 0 in a PM phase. We measure all the 3 × 3 rhomboids in the lattice to obtain the expectation value 〈string〉 along the path scanning μ at V = 0.9. The resultant data in Fig. 3d indeed reveal that inside the odd (even) Z 2 QSL phase, 〈string〉~−1(〈string〉~1), while inside the PM phase, 〈string〉 ≈ 0; the transitions are also seen to be first-order, in consistency with Fig. 3b, c.

The dynamical dimer spectra
One of the hallmarks of a QSL is its ability to support fractionalized excitations that cannot be created individually by any local operator. In this section, we focus on one class of such fractional excitations with magnetic Z 2 flux, i.e., the visons. Naturally, vison configurations with different fluxes will result in different dimer spectral signatures, thus realizing, in particular, the interesting phenomenon of translational symmetry fractionalization [10][11][12][13]37 , which can be further used to distinguish the PM and the even/odd Z 2 QSLs and make a possible connection to experiments. To this end, we compute the dimer spectra, obtained from stochastic analytic continuation of the Monte Carloaveraged dynamic dimer correlation function D(k, τ) with τ ∈ [0, β] (which can be viewed as the dynamical vison-pair correlation functions deep inside the Z 2 QSLs 37 ; more details can be found in the Supplementary Note 1). Figure 4a shows that in the odd Z 2 QSL phase, the gapped dimer (vison-pair) spectrum forms a continuum, and the dispersion minima are located at both the M and Γ points 35,37 . On the other hand, Fig. 4c illustrates that the dimer (vison-pair) spectrum deep inside the even Z 2 QSL is also a continuum but with minima only at Γ. These features are consistent with the expectation that the visons of the odd Z 2 QSL carry a fractional crystal momentum, whereas visons of the even QSL do not 12,37 . Note that for the single vison dispersion of an odd QSL, the locations of the minima are dependent on the chosen gauge 30,47,48 whereas the vison-pair spectrum is a gauge-invariant observable. For the even QSL, refs. 24,25 found that the minima of the mean-field vison dispersion occur at the three inequivalent M points in the Brillouin zone. Accordingly, one would then expect the vison-pair spectrum to exhibit a minimum at Γ (which is equivalent to 2M modulo a reciprocal lattice vector), in agreement with our numerical results. The arguments above apply generally to the dynamics of an odd/even QSL and should hold even at finite μ; similar behaviors have also been observed for the odd/even QSLs of the Balents-Fisher-Girvin (BFG) model 12 . In comparison, Fig. 4b presents the dimer spectrum inside the PM phase; here, there exists no clear continuum in the frequency  Fig. 1. All the data are simulated using β = L = 12. The upper-right labels in each panel represent the scaling factor for the intensities such that the five panels can be scaled onto the same color bar. In addition, the high-symmetry path for the spectra in Fig. 4 is also drawn in (a).
domain, indicating the lack of fractionalization of dimers into pairs of visons. Moreover, the overall dispersion is flat, which is consistent with the dispersionless S z spectrum in an S x -polarized state, such as in the transverse-field Ising model.

Discussion
In this work, we investigate a QDM with variable dimer density on the triangular lattice and uncover a plethora of interesting phases, including crystalline solids and two distinct classes of highly entangled QSL states hosting fractionalized excitations. Through detailed quantum Monte Carlo analyses, we explore the subtle interplay between these different phases and find the unique properties of their static and dynamic fingerprints.  ω ω ω Fig. 4 | Dynamical dimer spectra. The dimer spectra in the (a) odd Z 2 QSL in the limit of one dimer per site, corresponding to μ → − ∞ and V = 1 in Fig. 1, b PM phase with μ = 0, V = 0.9 and h = 0.4, and c even Z 2 QSL in the limit of two dimers per site, corresponding to μ → ∞ and V = 0.5 in Fig. 1. The dimer spectra exhibit continua in both (a) and (c), conveying the fractionalization of spins into visons. However, the dispersion minima in the two cases differ, being located at both M and Γ for (a) and only at Γ for (c), representing the translational symmetry fractionalization in the former and the lack thereof in the latter. In (b), however, the dimer spectrum is flat and displays less of a continuum in the frequency domain, consistent with a polarized PM phase. All the data are simulated at β = 200 on a L = 12 lattice, with the low-temperature T = 1/200 being necessary to overcome the small vison gap and the transverse field h. In particular, our results could find application to recent experiments with programmable quantum simulators based on highly tunable Rydberg-atom arrays, which have emerged as powerful platforms to study strongly correlated phases of matter and their dynamics. While our extended QDM differs from models of Rydberg atoms on the sites of the kagome lattice 28 in the precise form of the V interactions, the two systems bear resemblance in some of their phases. Specifically, the Rydberg model also displays the 1/6 staggered and nematic phases of Fig. 1, separated by a 'liquid' regime with no broken symmetry. These ordered phases can be mapped to the solid phases of a triangular-lattice QDM with either one or two dimers per site, which precisely constitutes our soft constraint. Appealing to the universality of phase transitions 28 , possible fates of the liquid state in the Rydberg model are then one or more of the phases obtained by interpolating between the 1/6 staggered and nematic phases in Fig. 1 for the present quantum dimer model: namely, the odd QSL, the PM, and the even QSL. These considerations highlight the potential utility of variabledensity dimer models in the experimental realm and provide a pathway to studying their rich physics.

Sweeping cluster algorithm
This is a quantum Monte Carlo method developed by the authors to solve the path integral of constrained quantum many-body models [37][38][39][40]49 . The key idea of the sweeping cluster algorithm is to sweep and update layer by layer along the imaginary-time direction, so that the local constraints (gauge fields) are recorded by update lines. In this way, all the samplings are performed in the restricted Hilbert space, i.e., the low-energy space. The original sweeping cluster QMC method 38,39 is designed for hard-constraint models, i.e., models in which the number of dimer(s) per site is fixed 37,50 . To solve our models in this work, we further improve upon the prior methods to be able to simulate a soft-constrained dimer model.
The Hamiltonian that we consider is given by Eq. (1) supplemented with the "soft" constraint that there can only be either one or two dimer(s) per site. The definition of winding numbers 39,51-54 for these two cases are explained in Supplementary Note 3.
Similar to the practice in Stochastic Series Expansion types of quantum Monte Carlo methods 55 , we separate the Hamiltonian into diagonal and off-diagonal parts. It is obvious that the t and V terms will not change the number of dimer(s) per site, but both the chemical potential μ and the transverse field term h would. Therefore, the Monte Carlo update will need to obey the soft constraint when we deal with the μ and h terms. We write the h off-diagonal term and the μ diagonal term as, where C is a constant to ensure that the corresponding matrix elements are positive. The label "d/o" indicates whether the operator is diagonal or off-diagonal, and l labels the links of the lattice. Although these two terms are single-link operators, they may break the soft constraint when considering neighbors, so we have to regard the single-link operator as a multi-link operator instead with all closest neighbors as shown in Fig. 5.
We can design the Monte Carlo algorithm to update vertices according to the soft constraint on the cells as shown in Supplementary Fig. 5. Since the original sweeping cluster method always obeys the constraints without changing the number of dimers per site, adding such considerations for the terms in Eq. (1) into the original sweeping cluster Monte Carlo method makes all samplings satisfy the soft constraint.

Stochastic analytic continuation
The main idea behind the stochastic analytic continuation (SAC) method 41,42,46,56 is to obtain the optimal solution of the inverse Laplace transform via sampling dependent on the importance of goodness. A set of imaginary-time correlation functions G(τ) can be obtained through the sweeping cluster QMC method first. The real-frequency spectral function and the imaginary-time correlation function are related by a Laplace transformation as GðτÞ = R 1 0 dωðe Àτω + e ÀðβÀτÞω ÞSðωÞ=π. We can inversely solve this equation by fitting a better spectral function. Assume the spectral function has a general form, S(ω) = ∑ i a i δ(ω − ω i ). We can obtain the optimal spectral function, i.e., the optimal choice of the set {a i , ω i } in the ansatz, numerically through sampling according to the importance of goodness of fit, with a simulated-annealing approach and with respect to the QMC errorbars of the imaginary-time correlation data G(τ). The reliability of such a QMC-SAC scheme has been extensively tested in various quantum many-body systems, such as the 1D Heisenberg chain 57 compared to the Bethe ansatz, the 2D Heisenberg model 44,58 in comparison to exact diagonalization, field theoretical analysis and neutron scattering spectra in real square-lattice quantum magnets, deconfined quantum critical points 58,59 and deconfined U(1) spin liquid phases with emergent photon excitations 60 , Z 2 quantum spin liquid models with fractionalized spectra 12,13,61 via anyon condensation theory, and the quantum Ising model with direct comparison to neutron scattering and NMR experiments 62,63 . We refer the readers to the technical descriptions available in the literature for detailed documentation of our QMC+SAC scheme.

Data availability
The data that support the findings of this study are available from the authors upon reasonable request. For the soft constraint of 1 or 2 dimer(s) per site, we have to consider all the neighbors when creating/annihilating a dimer on the central link. a When both the A and B sites have one dimer, one is allowed to create/annihilate a dimer on the center link. b It is forbidden to create a dimer on the center link when either the A or the B site already has two dimers. c It is forbidden to annihilate a dimer on the center link when either the A or the B site has only one dimer.